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Depending on the growth condition bacterial colonies can exhibit different morphologies. As 
argued by Ben- Jacob et al. there is biological and modeling evidence that a non-linear diffusion 
coefficient of the type D{b) = Dob k is a basic mechanism which underlies almost all of the patterns 
and which generates a long wavelength instability. We study a reaction-diffusion system with a 
non-linear diffusion coefficient and find that a unique planar traveling front solution exists whose 
velocity is uniquely determined by k and D = Do/D n where D„ is the diffusion coefficient of the 
nutrient. Due to the fact that the bacterial diffusion coefficient vanishes when b — > 0, in the front 
solution b vanishes in a singular way. As a result the standard linear stability analysis for fronts 
cannot be used. We introduce an extension of the stability analysis which can be applied to singular 
fronts, and use the method to perform a linear stability analysis of the planar bacteriological growth 
front. We show that a non-linear diffusion coefficient generates a long wavelength instability for 
k > and D < D c (k). We map out the region of stability in the D-fc-plane and determine the onset 
of stability which is given by D c (k). Both, for D — » and k — » oo the dynamics of the growth zone 
essentially reduces to that of a sharp interface problem which is reminiscent of a so-called one-sided 
growth problem where the growth velocity is proportional to the gradient of a diffusion field ahead of 
the interface. The moving boundary approximation that we derive in these limits is quite accurate 
but surprisingly does not become a proper asymptotic theory in the mathematical strict sense in the 
limit D — > 0, due to lack of full separation of scales on all dynamically relevant length scales. Our 
linear stability analysis and sharp interface formulation will also be applicable to other examples 
of interface formation due to nonlinear diffusion, like in porous media or in the problem of vortex 
motion in superconductors. 

PACS numbers: 5.40+j, 5.70.Ln, 6f.50.Cj 



I. INTRODUCTION 



A. Background of the Problem 



Recently the growth of bacterial colonies under differ- 
ent growth conditions has been the focus of attention 
of several groups in the physics community since it ex- 
hibits different elaborate branching patterns. For an ex- 
tensive review and entrance to the literature, see [1-3]. 
Already in 1989, Fujikawa and Matsushita [4] stressed 
that bacterial colonies could grow patterns similar to the 
type known from the study of physical systems such as 
diffusion-limited aggregation. A complete morphology 
diagram has been obtained for the colonics of Bacillus 
subtilis [1,5,6], where the important control parameter 
are the agar concentration which influences the diffu- 
sion of the bacteria as well as of the nutrient, and the 
initial nutrient concentration. It includes some interest- 
ing regimes such as diffusion limited aggregation, dense 
branching morphologies, Eden-like and ring patterns. Al- 
though the visual appearance of the patterns is very sim- 
ilar to those of physical systems, at the microscopic level 



their growth mechanism has to be different — the ques- 
tion then becomes whether or not these microscopic dif- 
ferences affect the overall large-scale pattern dynamics. 
For instance, the building units are bacteria which are 
themselves micro-organisms and thus living systems. To 
survive they have to cope with hostile environmental con- 
ditions which made them develop quite sophisticated co- 
operation mechanisms and communication skills, such as 
direct cell-cell interaction via extra-membrane polymers, 
collective production of extra-cellular "wetting" fluid for 
movement on hard surfaces, long-range chemical signal- 
ing, such as quorum sensing and chemotactic signaling, 
just to name a few. Different models have been proposed 
which include one or several of these mechanisms, and 
are able to reproduce the rich morphology diagram quite 
well. Instead of exploring the richness and diversity of 
the behavior of bacterial colonies, we want to concen- 
trate on the basic mechanism which underlies all these 
patterns. Since they appear as an interface separating a 
region occupied by the bacteria from a bacteria-free re- 
gion which propagates as the colony is expanding, we look 
for an interface model which includes a long wavelength 
instability. Although these models have been developed 
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and studied for pattern-forming, non-living systems such 
as crystal growth [7-9], where a sharp interface formu- 
lation is well justified, even at quite small length scales, 
here the existence of interface-type fronts is not obvi- 
ous from the start, but is something that should emerge 
from the continuum equations describing the dynamics.. 
Reaction-diffusion type models with a non-linear diffu- 
sion coefficient for the bacteria density have been argued 
to be a good candidate for being the proper starting point 
to analyze the instability mechanism since they were able 
to reproduce many aspects of the above mentioned mor- 
phology diagram [2,10-13]. 

The biological motivation that has been proposed for 
non-linear diffusion coefficient is the way bacteria move. 
Although there are different ways of moving we are in- 
terested in bacteria which swim by propelling themselves 
with their flagella in straight lines and change their di- 
rection in a random fashion by tumbling which can be 
described by a random walk. However, for the propelling 
mechanism to work a liquid with low viscosity is required. 
Since bacteria by themselves are able to secret this liq- 
uid, their presence is required to generate the lubricant 
layer necessary for diffusion. This behavior can be cap- 
tured qualitatively by a bacteria density dependent dif- 
fusion coefficient as has been proposed in particular by 
Ben- Jacob and co-workers in [11]. A consequence of it 
is that the branches of bacterial colonies are confined by 
a sharp envelope which is supported by the observation 
with optical microscopes [2,11]. 

However, we would like to not that the arguments sup- 
porting a non-linear diffusion model are still not conclu- 
sive and more of a qualitative nature. In addition it is 
clear that it does not appear to be relevant for the growth 
patterns at large agar concentrations where the bacteria 
are non-motile [6] , and the relevance for the regions where 
bacteria are motile is still under discussion. In this paper 
we will not address the question of the biological rele- 
vance of the model; instead we aim to contribute to the 
debate by working out the stability diagram and uncover- 
ing its essential dynamics. An additional non-biological 
contribution of our paper is that we introduce new meth- 
ods to deal mathematically with singular fronts. As we 
discuss below, this is likely to have implications in other 
subfields of physics. 

In passing, we also note that it has been shown re- 
cently [14] that if one extends the model by introducing 
an effective cutoff in the reaction term modelling the bac- 
terial growth, while keeping the bacterial diffusion term 
linear, one also recovers the type of front instability nec- 
essary to understand bacterial patterns. The motivation 
for such a cutoff would be simply the fact that bacte- 
ria are discrete entities, so that at some small density a 
continuum formulation breaks down. Both mechanism 
(nonlinear diffusion and discrete entity cutoff effects to 
continuum formulations) are not mutually exclusive and 
can be operative simultaneously, but the detailed studies 
of various models by a number of authors [2,10-13] sug- 
gests that the nonlinear diffusion mechanism is the most 



important one of the two [15]. 

We concentrate on the effect of a nonlinear diffusion 
coefficient here since in spite of the suggestion that a 
nonlinear diffusion coefficient is a possible mechanism to 
generate the complex morphology diagram, a clear un- 
derstanding of this instability mechanism is still missing. 
This is surprising since also from a mathematical point of 
view it is an interesting problem as it defines a new class 
of fronts which do show up in other systems with den- 
sity dependent diffusivity, such as porous media [16-18]. 
Furthermore, magnetic flux vortices in superconductors 
[19,20]. Clearly, understanding the similarities and dif- 
ferences between instabilities in magnetic flux patterns 
and the well-studied Mullins-Sekerka instability mecha- 
nism is clearly of importance. Considering the amount 
of work and attention there has been in the recent years 
to understand the mechanisms behind bacterial colony 
growth it might at first sight seem surprising that not 
even a stability analysis of planar fronts solutions has 
been performed. An important reason for this is that as 
the problem involves singularities: these make the stan- 
dard stability calculations break down, so new techniques 
have to be introduced to even perform the linear stabil- 
ity analysis. We have been able to resolve the problem 
and thus perform an explicit linear stability analysis of 
planar fronts which allows us to determine the regions of 
stability in parameter space. Our extension of the stan- 
dard stability calculation is not limited to the particular 
bacterial growth problem we focus on here. Instead it 
should be applicable to a large class of growth problems 
with singular fields, e.g., other problems which involve 
nonlinear diffusion, like the vortex patterns in supercon- 
ductors [19,20] just mentioned, should be amenable to 
the same type of analysis. 

In some limits, in particular in the limit that the bac- 
terial diffusion coefficient becomes much smaller than the 
one of the nutrient, the fronts in the models that have 
been studied become rather sharp. A second important 
question therefore is to what extent a moving boundary 
approximation, in which the front is viewed as a math- 
ematically sharp interface on the scale of the patterns, 
becomes appropriate — such approximations are often 
very helpful for analyzing pattern forming problems (see 
e.g. [21] for an application to dendritic growth and an 
entry into the vast "phase field model" literature) . Some 
steps in this direction for the bacterial growth problem 
were taken by Kitsunezaki [22] . We address this question 
in more detail in this paper and, quite remarkably, find 
that while in the limit of small bacterial diffusion a mov- 
ing boundary approximation is quite accurate it does not 
emerge as the lowest order description in a mathemati- 
cally well-defined limit. The reason for this is that even 
for small diffusion, the dynamically relevant length scales 
(i.e., those corresponding to unstable modes in the lin- 
ear stability calculation of planar fronts) are not all large 
in comparison with the front width. This result shows 
that bacterial growth problems with nonlinear diffusion 
of the type encountered in the porous media equation 
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[16,18] are mathematically in some crucial ways different 
from the standard type of growth problems. Physically 
their dynamics is closest to those of the so-called one- 
sided growth problems [7] . 



B. The Model 

Since we would like to concentrate on the basic mech- 
anism which generates a long-wavelength instability we 
confine our analysis to the most basic model of Ben- 
Jacob et al. [10,11], namely a two-dimensional reaction- 
diffusion system for the bacteria density b(r, t) with a 
nonlinear diffusion term, and the nutrient density n(r, t) 
with a linear diffusion term: 



db 

dt 
dn 

~dt 



= VD(6)V6 + /(n, b), 
= D n \7 2 n- g(n, b), 



(1) 
(2) 



with D n describing the diffusion constant of the nutrient, 
and 



D(b) = D Q b k 



(3) 



implying a bacteria density dependent diffusion coeffi- 
cient as was motivated before. For simplicity we assume 
the following reaction term 



f(n, b) = g(n, b) = nb, 



(4) 



which in chemical terms is like a bilinear auto-catalytic 
reaction: 



N + B ^2B. 



(5) 



Biologically it models that the bacteria B eat a nutrient 
N to duplicate themselves. This involves a conservation 
law and is clearly an oversimplification, since part of the 
energy is also used for movement and other metabolic 
activities. For the growth process we want to study here, 
this should not matter. For the same reason, we also leave 
out in this paper another biologically important feature, 
sporulation, a transition of motile bacteria into a station- 
ary state; this occurs if there is a deficiency of nutrient, 
which seem to play an important role in the later stage of 
the branching process. During sporulation bacteria stop 
normal activity such as movement and use all their inter- 
nal reserves to metamorphose from an active motile cell 
to a spore, a sedentary durable "seed" which is immotile 
and hence cannot participate to the diffusion process. 
The sporulation process can be included in the model by 
adding a term — nb on the right side of (1). Although the 
simulation by Kitsunezaki [22] indicate that this death 
term does affect the stability of planar fronts, we will not 
take it into account here since the most crucial ingredi- 
ent is the nonlinear diffusion coefficient of b as it assumes 
that without bacteria there is no diffusion. As we will sec 



this implies a front profile which goes abruptly to zero, 
with a divergent slope for k > 1. This characteristic is 
supported by experimental observations of some kinds of 
bacteria, where one observes a clearly defined envelope 
(such a comparison suggests a value of k of about one). 
The question we want to study now, is whether this kind 
of diffusion is enough to generate a long wavelength insta- 
bility. It should be noted here, that for k = the system 
has been studied by [23-25]. They showed that bilinear 
autocatalysis alone is not sufficient to destabilize a pla- 
nar front. Only in the presence of an autocatalysis term 
proportional to V 1 with 7 > 1 and D n > (3 c Dq where 
(3 C depends on the amount and order of autocatalysis a 
planar front is unstable toward long wavelength pertur- 
bations [15,26,27]. Thus, any instability we observe for 
k > is due to the nonlinearity in the diffusion term. 
By rescaling the diffusion constant D — Do/D n and re- 
placing f(n, b) and g(n, b) by (4), we obtain the following 
nonlinear reaction-diffusion system: 



db 
dt 
dn 
~dt 



D 



k + 1 

= V 2 n - nb, 



V 2 b k+1 + nb, 



(6) 
(7) 



which contains two parameters, D the rescaled diffusion 
constant, and k describing the nonlinearity and the stiff- 
ness of the front — in writing the above equations, we 
have used the freedom to choose appropriate time and 
length scales, and to rescale the fields n and b appropri- 
ately to set all other prefactors equal to one. We will 
be interested in front solutions of this equation where far 
ahead of the front the nutrient field n — > 1; as we will 
discuss in more detail below, this asymptotic value is 
also immaterial, as the problem with another asymptotic 
value can be rescaled to our problem with a rcnormalizcd 
value of D. 

A nonlinear diffusion behavior like in (6) also arises in 
the so-called porous media equation [16-18]. There is a 
vast literature on this equation [16,18]; for us, the essen- 
tial feature is that it gives rise to moving front solutions 
with compact support, i.e., for which the field b is zero in 
some regions of space. At the point where b vanishes, it 
does so in a singular way, and this invalidates the usual 
linear stability analysis. 



C. Overview of methods and results 

For the reader not interested in the mathematical de- 
tails of the derivation, we now summarize the main re- 
sults of the analysis. The model (6)- (7) has two homoge- 
neous states: a stable solution (q,, 0) in which only bac- 
teria are present, and an unstable solution (0, c„) with 
only nutrient. Thus, we can study the propagation of 
the stable state (q,,0) into the unstable one (0,c„), im- 
plying for our system the propagation of the bacteria field 
into the nutrient field. To study such a propagation we 
look for one-dimensional traveling front solutions which 
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appear for a system with initial conditions in which the 
system is in the unstable state and a small perturbation 
at x — > — oo starts to invade it. Assuming that the front 
propagates with a steady velocity v, we can reformulate 
the model in a co-moving frame which reduces (6)- (7) to 
a one dimensional system of ODE's which is much easier 
to analyze. Its solution will be found numerically by a 
shooting method as will be explained in section II. 

We find that there generally is a clearly defined unique 
reaction front, of which b vanishes with a diverging slope 
for k > 1 (see Figs. 5 and 6 below). The qualitative 
features of these fronts are consistent with the earlier 
simulation results of Ben- Jacob et al. [10,11] and can be 
traced back to the nonlinear diffusion. The characteris- 
tic singular behavior of the front makes the study of the 
problem mathematically and numerically challenging and 
intriguing. The solution provides us with a unique veloc- 
ity, which depends on D and k and is shown in Fig. 1. 
More detailed plots of the behavior of the velocity as a 
function of D and k are presented later in section II of 
this paper. 








FIG. 1. The front velocity as a function of D and k, as 
determined from the analysis in section II. 




FIG. 2. Perturbed front profiles of bacteria densities. The 
front propagates into the :r-direction, and has a sinusoidal 
modulation in the y-direction. 

In order to study the stability of the front which is the 
content of section III we have to perturb the planar front. 
Due to the singular behavior of the planar front a per- 
turbation of the front is not only a simple perturbation 
in the fields b and n but also in the geometry of the front 
as sketched in Fig. 2. 

Our stability analysis implements the idea that a 



proper Ansatz consists of two contributions, a perturba- 
tion in the line of the singular front, and the perturbation 
in the fields away from the singular line. Both these con- 
tributions have to be determined sclf-consistently. For 
k > we observe that for D < D c (k) the planar front 
is unstable and has a long wavelength instability. Thus, 
a nonlinear diffusion coefficient together with a bilinear 
autocatalysis-type reaction term are sufficient to gener- 
ate a long wavelength instability. For D > D c (k) the pla- 
nar front is linearly stable. Hence, in the D-fc-parameter 
space there exist regions of stability and instability of a 
planar front. We determine these regions in two different 
ways, one by performing numerically a linear stability 
analysis (LSA) as is done in section III. A, the other by 
an expansion for small growthrate u> and wavenumbcr 
q around the planar profile as is done in section III.C. 
Fig. 3 shows the stability diagram as a function of D and 
k. 
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FIG. 3. Stability diagram for parameters D and k. Filled 
circles show where the region of stability of planar fronts starts 
as determined by a numerical linear stability analysis, filled 
diamonds show the same boundary as obtained from the solv- 
ablity formula for d 2 u}/dq 2 \ q= o derived in section III.C. For 
k = 0.5 both methods give the same value up to the size 
of the symbol. The solid line is there to guide the eye, the 
dashed line hints at the fact that while we expect the line of 
D c (k) to approach the origin we do not know the precise an- 
alytic behavior of D c (k) for k — > 0, since for k — the planar 
front is stable for all D [24,23]. The two crosses represent 
the simulation performed by Kitsunezaki [22]. For D = 0.2 
the front in these simulations was unstable which is consistent 
with our analysis. For D = 1.0 the planar front was stable, 
which does not agree with our analysis. The probable cause 
of this apparent discrepancy is discussed in the main text. 

Filled circles show the onset of the region of stability 
of planar fronts as determined by a numerical linear sta- 
bility analysis, filled diamonds show the same boundary 
as obtained from the exact expression for d?u}/dq 2 \ q= o 
derived in section III.C. Both methods give results which 
are in very good agreement with each other, as they 
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should. The solid line is there to guide the eye, the 
dashed line hints to the fact that while we expect the 
line of D c (k) to approach the origin we do not know the 
precise analytic behavior of D c (k) for k — > 0, since for 
k = the planar front is stable for all D [24,23]. We 
will not analyse the precise behavior in the limit k — > 
in detail, both because it does not appear to be of prac- 
tical relevance, and because the model is very sensitive 
to slight changes in this limit: an effective cut-off which 
arises for discrete particle effects turns the model weakly 
unstable [14], but a continuum model with a different 
reaction term has the same effect. In particular, if we 
change the reaction term nb in (6), (7) to nb 1 , then for 
any 7 > 1 we expect for the limit k — > D c to be finite; 
in other words, for 7 > 1 the stability boundary crosses 
the D-axis at a nonzero value of D. For 7 = 2, it is in 
fact known that D c (k = 0) « 0.34 [24]. 

The two crosses in Fig. 3 represent the simulation per- 
formed by Kitsunezaki [22]. Whereas for D = 0.2 his 
planar front was unstable which is consistent with our 
analysis, his planar front for D = 1.0 appeared to be sta- 
ble, in apparent contradiction with our results. However, 



the simulations were done for a system of width L y 
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and up to time t = 200. From our results for the disper- 
sion relation for k = 1 and D = 1 which is very similar to 
the one shown in Fig. 10 in section III, we find that the 
characteristic length scale of the fastest growing mode 
is L m « 31, while the associated characteristic time for 
this fastest growing mode is approximately t m = 520. 
Hence, it is likely that the system width is too small and 
the simulation time too short to observe the instability. 
It would therefore be useful to redo the simulation for a 
bigger system and longer times. In fact, this illustrates 
the difficulty of using simulations alone to study the sys- 
tems, especially if only a few parameter values can be 
studied over a limited time range and system size. On 
the other hand, our explicit stability analysis allows us 
to map out the phase diagram in a relatively straightfor- 
ward way. 

In section IV we map the system with a moving bound- 
ary approximation to a sharp interface problem guided 
by the success this approach had in analyzing and under- 
standing the Mullins-Sckerka instability mechanism [7], 
the long wavelength instability associated very generally 
with diffusion-limited or Laplacian growth processes. We 
obtain by a multiscale expansion equations for b and n 
which arc valid in the outer bulk fields, and which are 
connected by boundary conditions. The boundary con- 
ditions are obtained by using solvability type arguments 
to integrate out the internal degrees of freedom of the 
inner reaction region. As was already mentioned before, 
the moving boundary approximation is closest to the so- 
called one-sided growth models and is quite accurate for 
small D, but it never becomes mathematically correct 
in the limit D — > for all dynamically relevant length 
scales. 



II. PLANAR FRONT 

There exist two trivial homogeneous solutions: 
n(x,t) = c n ,b(x,t) = 0, which implies some constant 
food level and no bacteria. This state is unstable since 
any amount of bacteria will be enough to let the bacte- 
ria density grow. The other trivial homogeneous state is 
n{x,t) = 0,b(x,t) = Cb, which assumes a constant bac- 
teria density and no food. This state is stable in the 
present model without sporulation. In addition there ex- 
ist a steady-state solution in which the stable state (q,, 0) 
propagates with a constant velocity v into the unstable 
state (0,c„), implying the propagation of the bacteria 
field into the nutrient field. Starting from an initial con- 
dition in which the unstable nutrient state is perturbed 
by a small amount of bacteria at the left, the bacteria 
field invades the nutrient state in the form of a well de- 
fined reaction front propagating to the right. Since we are 
first interested in a planar front, we can restrict ourselves 
to one dimension. To obtain the uniformly translating 
front solution it is convenient to express the reaction- 
diffusion system in a co-moving frame in which the new 
coordinate £ travels with the velocity i>o of the front, 
£ = x — vot. The temporal derivative then transforms 
as dt\ x — dt\t — vod^\t- For a front translating with uni- 
form velocity i>o, the explicit time derivative vanishes and 
(6)-(7) reduces to: 



D 



d 2 b k+1 db , „ 



d 2 n dn 

de +Vo d~r nb 



0. 



(8) 
(9) 



This is a system of two ODE's of second order. The 
boundary conditions at £ — > ±00 are given by the two 
homogeneous states. By choosing a right-moving front 
we obtain as boundary conditions at £ — > — 00 the stable 
state: 



K£ -> -°°) = c <" d H b (t, 
n(£ — > — 00) = 0, d^n(^ 



-00) =0, (10) 
-00) = 0, (11) 



which invades the unstable state given at £ — > 00: 

6(£ -> 00) = 0, c%&(£ -> 00) = 0, (12) 
n(£ — > 00) = c„, d|7i(£ — > 00) = 0. (13) 

As mentioned before, the system simplifies extremely in 
the region where &(£) = 0. By choosing the origin £ = 
in such a way that for positive £ &(£) = 0, the system 
(8)-(9) reduces in the positive ^-region to: 



KO = 0, 

d 2 n dn 

de +vo d-r°- 



(14) 
(15) 



which is a linear ODE for n which can be solved analyt- 
ically and is given by: 
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n(0 = c n - c exp(-w o O. 



(16) 



where cq > is determined by the full problem. Hence, 
the system can be divided into two regimes, the first be- 
ing £ > given by (14)-(15) which can be solved ana- 
lytically, and the second being £ < which contains the 
full nonlinearity. Both regimes are connected via their 
common boundary condition at £ = 0. Hence, it is suf- 
ficient to study (8)-(9) for £ < 0, for which we still have 
to determine the behavior at £ — > which we will obtain 
by studying the local behavior of the bacteria density b 
and the nutrient density n as £ approaches zero from the 
left. Since the bacteria density b is a physical quantity, 
we assume it to be continuous. Moreover, (9) then im- 
plies that n and its derivative at the boundary have to be 
continous as well. Hence, we obtain for n the boundary 
condition at £ = 0: 



n(0) 
dn 



CO, : 



= V C . 



(17) 
(18) 



In the introduction we have already discussed that the 
bacteria density b shows a singular behavior for £ — > 0. 
This is due to the fact that the prefactor of the highest 
derivative in the 6-equation contains a factor b k , which 
vanishes as b — > 0. This allows b to become singular near 
£ = 0. As is well known (see e.g. [28]) at such a regular 
singular point one expects a behavior for b of the type 
[29]: 



(19) 



Substituting this Ansatz into (8) we obtain: 

Da[a(k + 1) - 1] A k+1 (-^) a (*+ 1 )- 2 (20) 

-voAai-O"' 1 -n a A(-O a =0. 

To fulfill this equation the dominant terms in £, which 
are the first and second terms, have to cancel. This de- 
termines A and a to be: 



A 



kvp 
D 



l/k 



a= k- 



Hence the bacteria density profile vanishes as 



for £ -> 0, 



(21) 
(22) 

(23) 



which implies that its derivative db/d£ diverges for k > 1 
as 

^--^-O 1 '*" 1 for ^0. (24) 

Hence, we are left to study (8)-(9) for £ < with the 
boundary conditions (10), (11), (17), (23) and (24). 



Due to the fact that we chose /(n, 6) = g(n,b), a con- 
servation law is underlying the system (6)- (7), express- 
ing that all food is transformed into bacteria, i.e. that 
Cb = c n . The conservation law allows us to reduce the 
order of our system of ODE's by one. Hence, by adding 
(8) and (9) and integrating over space from — oo to £ we 
obtain: 



D d 2 b k+1 
k + 1 d£ 2 



db 

vo-r: + nb 
d£ 



0. 



dn D db k+1 



(25) 
(26) 



Note that (26) immediately implies c& = c„ since the 
derivatives all vanish at £ ± oo. This just expresses that 
food is converted into bacteria in this simplified model. 

The one-dimensional front profile governed by (8)-(9) 
can be represented by a heteroclinic orbit in the (6, d^b, n) 
phase space connecting the two steady states correspond- 
ing to the boundary conditions (10) to (13). Due to the 
possibility of solving the system of ODE's analytically in 
the positive £ region the front profile can be found by ap- 
plying a standard shooting method to the region £ < 0. 
By shooting from £ — > — oo along the unstable manifold 
and requiring it to connect to the trajectory flowing into 
the singular origin with the boundary conditions (17), 
(23) and (24) a relationship between the velocity v and 
the boundary condition n(£ = 0) is uniquely selected. 
The existence of a unique front solution is also consistent 
with so-called counting arguments for the dimensions of 
the stable and unstable manifolds of the fixed points of 
the flow. On the left, for £ — > — oo, there is only one un- 
stable mode leaving the homogeneous fixed point, which 
then fixes n and d^n at £ — > completely. Matching at 
£ = to the positive £ solution for n can only be done on 
a line in the n — c^n-plane, since the n-solution is an ex- 
ponential. Hence, changing vo so as to match both fixes 
vq completely. 

As we already anticipated at the end of section LB, we 
henceforth choose c n = 1, and hence c& = 1: By appro- 
priately rescaling £, vo and the n and b fields, any other 
choice for c„ can be transformed to the case with c n = 1 
with renormalized diffusion coefficent Dr — Dc k . The 
uniquely determined front velocity is hence essentially a 
function of D and k only 



v = v (D, k). 



(27) 



We shall now study the behavior of the front profiles 
and of vq{D, k) in more detail by a combination of ob- 
servations from the numerical calculations and of simple 
analytical arguments. Many of these arguments can eas- 
ily be formalized by asymptotic analysis or by reducing 
the equations in certain limits to simpler ones, but we 
shall refrain from doing so explicitly. 

Fig. 1 gives an idea of the functional dependence of the 
velocity v on D and k. Fig. 4 displays that for small D 
the velocity is linear in D: 
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v « a(k)D (D < 1), (28) 

where a(k) is a proportionality constant which decreases 
with increasing k. 




FIG. 4. Dependence of the planar velocity vo on D for 
k = 1. The inset shows that for _D — ► the velocity ap- 
proaches zero linearly. 

This proportionality of v a with D for small D is simply 
a consequence of the fact that the propagation of the pro- 
file for small 6 is governed by the balance of the nonlinear 
diffusion with the v^db/d^ term. 




FIG. 5. Bacteria and nutrient density profiles for different 
D and fixed k — 2. 

Fig. 5 shows the dependence of the profile on D for 
k = 2. With decreasing D the interfacial thickness W 
decreases, whereas the diffusion length of the nutrient 
density l n increases since 



in = 1/vo 



(29) 



as seen from (16). Hence, with decreasing D there is a 
separation of scales between the diffusion length i n and 
the interface width. 




FIG. 6. Bacteria and nutrient density profiles for different 
k and fixed D = 0.3. 



Fig. 6 shows the dependence of the profile on k for fixed 
D, here D = 0.3. It demonstrates that with increasing k 
the interfacial region decreases and sharpens. 

At first sight, both Figs. 5 and 6 suggest that for small 
D or large k a moving boundary approximation might be- 
come appropriate. However, the behavior is rather sub- 
tle, and to prepare for a full discussion of this issue in 
section IV, we analyze the scaling of the front profiles in 
some more detail. 



0.35 



0.25 



0.5 



1.5 



FIG. 7. Interfacial thickness as a function of D for fixed 
k = 2. W is the distance from the origin to the point at 
which b = 0.5. 
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shown). Furthermore, for vo small and b w 1, Eq. (9) for 
n reduces to d 2 n/d£ 2 — n = 0, showing that n(£) decays 
to the left as e~^. In other words, n decays into the bac- 
terial zone on a length scale of order unity. Through the 
coupling term in Eq. (8), this also means that b decays 
to 1 towards the left on a scale of order unity — this is 
actually visible in Fig. 6. Thus, even though for large k 
b rises to values close to 1 on exponentially small length 
scales W, the scales over which b and n decay to their 
asymptotic values are actually of order unity. 



III. LINEAR STABILITY ANALYSIS OF 
PLANAR FRONTS 



FIG. 8. Interfacial position as a function of k and for fixed 
D = 0.3. W is the distance from the origin to the point at 
which b — 0.5. 

To quantify the behavior of the interfacial thickness 
W as a function of k and D let us first measure the 
thickness at which the bacteria density reaches the level 
b{W) = bw = 0.5. Fig 7 shows how the interfacial thick- 
ness approaches a finite thickness as D approaches zero, 
and Fig. 8 how W approaches zero with increasing k. 
Both dependencies can be understood by inverting (23): 



W 



v k bw - 



(30) 



Since v is proportional to D for small D, the interfacial 
thickness approaches a constant W : 



W^W = 



1 



a{k)k 



h k 



w 



(31) 



which depends only on k and the chosen interfacial value 
bw- With increasing k, W$ decreases, and vanishes for 
k — > oo; indeed, for not too large values bw, we have 



b k « cxp {-k\ ln&vH), 



(32) 



so that Wo vanishes exponentially Note finally that 
Fig. 8 indicates that W becomes large as k — > 0; this 
indicates that the behavior of the model for k <C 1 is 
quite different from that in the regime k of order 1 or 
larger, on which we will concentrate. 

So far, we analyzed the width between the point where 
b reaches some fixed value bw < 1 and the point where 
b vanishes. In the limit D — > this width remains finite, 
while for k — > oo the width measured this way vanishes. 
However, for addressing the question whether a sharp in- 
terface formulation can capture the essential behavior, it 
is also important to analyze how b approaches the asymp- 
totic value 1 for large k. When k is large, we see that n(£) 
becomes small in the interfacial zone. In fact, it is easy 
to convince oneself that the self-consistent scaling behav- 
ior of Eqs. (8)-(9) for £ < is n(£) - 1/fc, v ~ 1/k for 
large k, and this is born out by our numerical results (not 



A. Dispersion relation 

To study the linear stability of the planar front, wc 
have to perturb the front. Due to the singular behavior 
of the planar front the dynamically relevant perturba- 
tions are not just simply perturbations in the fields b 
and n but also in the shape of the singular line where 
6^0. Since we only study the linear stability, we al- 
low the perturbations to be complex and wc can focus 
on a single mode with wavenumber q and amplitude e by 
writing 

h(y,t) = eexp(iqy + wt), 

We take this function h to be the modulation of the po- 
sition of the line where the bacterial front vanishes, as 
indicated in Fig. 2. To be concrete, we now write b and 
n as 

6(6 tf,t) = 6o(£ + %.*)) + 

e6i(£ + %.*))<sxp(*<W + wt), (33) 

n(Z,y,t)=n {Z + h{y,t)) + 

em(^ + h(y,t))exp(iqy + ut). (34) 

where (bo, no) is the planar front solution determined in 
the previous section. This ansatz is the crucial ingredient 
that makes our stability analysis possible. The standard 
perturbation approach would amount to writing the per- 
turbed field b as b = b (£)+eb 1 (t;)e iqy+ " t ; such an Ansatz 
works only if &o(£) is smooth enough that its derivative 
remains finite — here, because of the singular behavior of 
bo, this standard approach fails. We therefore shift both 
the position of the singularity line of bo and of b\ , where 
&i and n\ are the corrections to the bacterial profile and 
nutrition field as a result of this modulation. In order 
that perturbations are arbitrarily small as e — > so that 
we can linearize the equations, we clearly need to have 



— bounded, ^ bounded. 
n &i 



(35) 



Moreover, of course, b\ and m should be continous twice 
diffcrentiable functions away from the singular line. 
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For the analysis, it will be convenient to introduce the 
locally co-moving frame 

C = x - v t + h(y, t) = £ + h(y, t) 

in terms of which the fields can be written as: 

6(6 y,*) = 60(C) + e&i(C)ex P (ujy + wi), (36) 
n{(, y, t) = n (C) + eni(C) CX P (ky + wt). (37) 

Upon linearization of the dynamical equations (6)-(7) 
about the uniformly translating solution (&o(0> n o(0)j 
we then get 



61 

m 



k+iJ y 








61 

m 



dbg 



(38) 



where / = 6q +1 and where the prime refers to a differ- 
entiation with respect to 6 . The terms proportional to 
dbo/dQ and dno/d( on the right result from the modu- 
lation h of the singular line about the line £ = in the 
argument C of 60 and no- Finally, the operator C is given 

by 



£11 



p d 2 

k + ldC 2 
k + 1 1 d( 2 



(/'•) 



vo 



n 



£12 

C21 = 

£■22 = 



D 

fe + 1 

60, 

-no, 

d 2 



f" 



dbo 



d 

D f „db Q \d 

ITi f -dc +Va )d( + 
d fl ,d 2 b Q 



(39) 



k + V dC, 2 



d ^+v --b . 



(40) 
(41) 

(42) 



Note that the eigenvalue equation (38) is an ODE prob- 
lem in terms of the variable £, in the same way as it is in 
the standard linear stability calculations [30] . 

Let us pause for a moment to reflect on the difference 
with the usual stability approach a bit more. Since the 
translational mode (e\6 , c^n ) is the right zero eigen- 
modc of £, 



This is precisely the linear equation one gets if one 
starts with the usual linear stability Ansatz b = 60(C) + 
61(C) exp(i<?y + tot), n = n (£) + ni(£) exp(iqy + cot) in 
terms of £ rather than ( as the variable. While at this 
level the two problems appear to be the same, their inter- 
pretation is not. When we write the perturbed problem 
in terms of the shifted coordinate ( and require 61/60 to 
remain bounded, then clearly (44) shows that the variable 
61 is more singular than 60 — in particular the singular 
behavior of 61 is that of db /d(. In other words, 61/60 is 
not a small perturbation, instead it diverges! Of course 
it is simply due to the fact that one can not represent 
a shift of the singular line with a small perturbation in 
terms of fields which vanish at £ = 0. The Ansatz we 
make in terms of the variable C, on the other hand, docs 
represent a proper shift of this line; it can be thought of 
as a suitable resummation to capture this. 

Let us return to the problem of solving for 61 (£) and 
ni((). Again we can split up the problem into two sepa- 
rate regions, since for ( > the problem simplifies to: 



61 =0, 



(46) 



d 2 n x 
d( 2 



v^-^ + q^n^iLO + q 2 )^. (47) 



dni 



which is a linear ODE in m which can be solved analyt- 
ically: 

m = (c n - c )v exp (-vqQ + d cxp (-AC), (48) 



with A = (v — \Jv 2 + 4(w + q 2 ))/2 and with d a some 
constant which is undetermined at this stage (co and c n 
are parameters of the solution (16) for no). This solution 
is connected to the negative C region via the boundary 
condition at ( = which determines do- To obtain the 
boundary condition at ( = 0, we analyse the local be- 
havior of 61 and n\ as ( — > from the left. Since ri\ and 
its derivative are continuous across ( = 0, m and c^m 
obey at C = 0: 



m = (c„ - c )u + do, 
d c ni = -(c„ - c )uo - do A. 



(49) 
(50) 



In view of our requirement (35) that 61/60 remains 
bounded, it is natural to assume that 61 vanishes as 



we see that if we introduce 
dbo 



61 = 61 + 



ni = m + 



dn 



d( ' ' ' ' d( ' 
in (38), then these new variables obey simply 

h \ _ (00 + ^f'q 2 \ fh 




uo + q z 



(43) 



(44) 



(45) 



61 = B(-CY 



(51) 



Indeed, by inserting it into (38) we straightforwardly ob- 
tain from the asymptotic behavior (19) for 60 



B = 



loA 
kv 



to ( kvo 
kvo~ [~D 



l/k 



<4 

Hence, for ( — > 0: 



(52) 
(53) 
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MO = -£<-<>'" 



kv 



(54) 



so that 



bo(0 



kv 



for £ -> 0, 



verifying that 61/60 remains finite. Hence a solution 
61 which vanishes according to (51) does obey the re- 
quirement that perturbations are small everywhere. The 
boundary conditions at ( — > — 00 are given by: 



61(C) -o, 9^i(0^o, 

ni(C)->0, a c ni(C)-0, 



(55) 
(56) 



since all perturbation should vanish at ( — > — 00. 

The linear dispersion relation is obtained by solving 
(38) for different q with the shooting method. By shoot- 
ing from ( — > — 00 along the unstable manifold and 
matching it to the trajectory leaving the origin with 
the boundary conditions (49), (50) and (54) we obtain 
a unique u as a function of D, k and q. At the same time 
do is determined. Counting arguments for the multiplic- 
ity again support the uniqueness of w. A numerical dis- 
persion relation was obtained for k = 0.2,0.3,0.5, 1,2,3 
and 5 and different D. For a fixed k the dependence on 
D of the dispersion relation is qualitatively the same for 
all k. Fig 9 shows the dispersion relation for k = 2 and 
different D. 

There is a long wavelength instability for all D < 
D c (k), whereas all modes are stable for D > D c (k). As 
D decreases below D c the growth rate of the unstable 
modes starts to increase as does the range of wavenum- 
bers which are unstable. At the same time both q m , the 
wave number which corresponds to the maximum growth 
rate, as well as q c , the wave number for which u — 0, shift 
with decreasing D to larger wave numbers. By decreas- 
ing D even further we observe that the growth rate starts 
to decrease again which is due to the fact that the whole 
dynamics of the front is slowing down as we decrease D. 
Note, however, that as D becomes small, the range of 
unstable wavenumbers does not vary appreciable: q c is 
roughly constant. 
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FIG. 9. Dispersion relation for k = 2 and different D. For 
D < D c the planar front is unstable for q < q c whereas for 
D > D c it is linearly stable for all q. 



0.02 



0.01 



CO 



-0.01 







• k = 0.2 

• k = 0.3 
k = 0.5 

k = 2 

k = 3 

k = 5 


■ .... ' x 








\ " N ' s 

. • \ 


\ \ 

\ 

\ 



0.2 



FIG. 10. Dispersion relation for D = 0.3 and different k. 

The dependence of the dispersion relation on k is 
shown in Fig. 10. We know that for k = the planar 
front is stable. For small k the front starts to be unsta- 
ble for long wavelength perturbations. With increasing k 
the range of unstable modes is increasing as is the growth 
rate. However at k > 1 the growth rate starts to decrease 
again which is again due to the fact that the whole dy- 
namics of the front slows down as k is increasing. q m 
shows qualitatively the same behavior as q c . It starts to 
shift with increasing k to shorter wave length, stays con- 
stant for k — 0.5 to k — 3 and then decreases again to 
longer wave length. Fig 3 shows how D c depends on k. 
With increasing k, the transition value D c increases, thus 
implying that with increasing k the region of instability 
is larger. For large k the value of D c (k) appears to be 
linear in k. 

One general noteworthy feature of our results is that 
the growth rate of the most unstable mode as well as the 
corresponding wavenumber q m are generally rather small. 
As we discussed already in section I, this may the reason 
that Kitsunezaki [22] appears to observe a planar stable 
interface in the region of the phase diagram where planar 
interfaces are unstable according to our calculation. 



B. Comparison with the Mullins-Sekerka instability 

The dispersion relation of the planar bacterial fronts is, 
for D < D c (k) and away from the instability line D c (k), 
very similar to the so-called Mullins-Sekerka dispersion 
relation 



u M s = v \q\(l - d e t hq 2 ) 



(57) 



that one derives for perturbations of a planar crystal- 
lization interface [7]. In this case, £ t h = D t h/vo is the 
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thermal diffusion length (the analogue of our nutrient dif- 
fusion length £ n ), and do is a microscopic surface-tension- 
like length which measures the strength of the curvature 
corrections to the interface. We shall see later in section 
IV why this analogy is justified, but it already shows us 
here something interesting: As D — * 0, Vq vanishes pro- 
portional to D. In this limit £th diverges just like t n does. 
Hence, from the observation that the range of unstable 
modes remains finite in this limit, and hence that the 
term analogous to d^lth remans finite, we can immedi- 
ately conclude that the "effective surface tension" of our 
bacterial fronts, the analogue of do in (57), should scale 
as D for small D. 
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FIG. 11. Sketch of a perturbed front propagating from the 
left to the right. The arrows drawn with a full line indicate 
the diffusion flow of nutrient, on the front side of the inter- 
face, those drawn with a dashed line the diffusion current of 
bacteria. At a protrusion into the nutrient region, the nutri- 
ent diffusion is enhanced while the bacterial diffusion current 
is suppressed. There are hence two competing effects, whose 
relative strength depends on D. 

That a propagating, planar reaction-diffusion front 
shows a long wavelength instability for small D but is 
linearly stable for all q for D > D c , has been observed 
and explained before (see e.g. [24]), and can be under- 
stood in the following way. Let us consider a perturbed 
front moving to the right as sketched in Fig. 11. At 
a protrusion into the nutrient side of the interface, the 
nutrient gradients are compressed and hence the nutri- 
ent diffusion is enhanced. The "feeding" of the interface 
from the nutrient side is hence enhanced there, and this 
tends to make such protrusions grow larger in time. On 
the other hand, as the dashed arrows indicate, the bac- 
terial diffusion flow from the back side towards the in- 
terface is reduced at such a protrusion — this tends to 
reduce the growth of such protrusions, and hence to sta- 
bilize the interfacial perturbation. The relative strength 
of the two effects is determined by D, the effective diffu- 
sion constant of the bacteria behind the interface. When 
D > D c (k), the stabilizing effect from the back side wins, 
for D < D c (k) the destabilizing effect on the front side 
dominates. Even the effect on k can be understood in 



this context. The effective diffusion coefficient is given 
by D — Db k , which lowers the effective diffusion coeffi- 
cient in the interfacial region where b < 1. Hence, the 
bigger k the smaller the effective diffusion constant in 
the interfacial region, the longer the destabilizing effect 
of the nutrient can prevail. When k decreases towards 
zero, the stabilizing bacterial diffusion extends more and 
more towards the front region [31]. 

As we pointed out above , in the limit D <C D c (k) the 
instability is very much like the classical Mullins-Sekerka 
instability of a crystal-melt interface. As D increases to- 
wards D c this connection breaks down because the sta- 
bilizing diffusion from the back-side becomes important 
within the interfacial zone: There is then no clear sep- 
aration anymore between an interface and the regions 
before and behind the front [see also section IV. C for 
further discussion of the behavior for D near D c (k)]. 

Of course, the competition between the stabilizing ef- 
fect of the diffusion gradient on the back side and the 
destabilizing effect of the gradient on the front side of 
the interface shows up in crystal growth during transient 
regimes and can be understood along the lines of the 
Mullins-Sekerka stability analysis [7]. A most amusing 
and dramatic illustration of this was observed recently in 
experiments on the melting of polarized 3 He [32]; there 
the instability sets in only after a very long transient be- 
cause the diffusion cocfficnt on the back side is very much 
bigger than on the front side; as a result, as long as there 
is a transient gradient on the back side, the melting in- 
terface remains stable. 



C. Onset of instability 

As we found above that the instability that occurs 
when D decreases below D c {k) is a long-wavelength q = 
instability the critical line D = D c (k) is the line where 
d(jj/d(q) 2 \ q= o = 0: to the right of this line in Fig. 3 this 
derivative is negative and to the left of it it is positive, 
so that lo > for small q. Since the translational mode 
q = is the eigenmode of C with eigenvalue u> — 0, we can 
investigate the behavior of the tu-q 2 curve in the vicinity 
of the origin by the following expansion [25]. Because the 
q = mode is a translation mode with zero eigenvalue, 
u) is small and of order q 2 when q is small. Moreover, 
&i = and ni = for q = 0, and so for small q, b\ and 
ni are both of order q 2 too. In (38) this implies that for q 
small, the terms on the right hand side involving b\ and 
ni are of order q 4 . To order q 2 , we therefore get 



bi 

"i 












(58) 



which is exact to order q 2 . Since C has a zero eigenvalue, 
we can apply the solvability condition by requiring that 
the inner product with the left zero mode of L vanishes: 



11 



*1 

*2 



k+lJ V 




<J \ I Fit 



u) + q 2 




(59) 



Here and ^>2 are the components of the left zero mode, 
i.e., of the right zero eigenvector of the adjoint matrix op- 
erator £*: 



C* 



D ci 
k+lJ 



l) 2 



bo 



a 2 



-no 

-v -§£-b Q 



(60) 



Upon rewriting (59) as 



f°° f T db Q dn 



D db dn 

*1 f + * 2 

k + i 1J d( 2 d( 



(61) 



and taking the limit q 2 — > 0, this leads to the required 
exact relation for the onset of the lateral instability: 



dw 



d(q 2 ) 



9=0 



dbo , \u 

1 ac ^ 



dn 

2-ac 



f) 



(62) 



Planar fronts change stability when the integral in the 
numerator of (62) changes sign. 

Since C is non-hermitean, there is no obvious relation- 
ship between the zero right eigenmode of C and its ad- 
joint C* . To find the zero right eigenvector of the adjoint 
operator C* we have to impose appropriate boundary 
conditions on the left eigenmodes too. Generally, the 
boundary conditions of the functions in the left adjoint 
space are obtained from the definition of the adjoint op- 
erator, in that for all functions $ we consider we need to 
have 



/oo />oo 
dC*(£$) = / dC(£*#)$. 
-oo J — oo 



(63) 



In general, when the partial integrations are done so as 
to obtain C* from C, we obtain boundary terms; the re- 
quirement that these vanish give the boundary conditions 
on the adjoint functions In the present case, since the 
functions on which our operators are working arc defined 
on the infinite interval (— oo, oo) for m and its related left 
component ^2, and on the semi-infinite interval (— 00, 0] 
for b\ and ^1, we find that the appropriate boundary 
condition for the adjoint functions is that ^2 should 
stay bounded as ±00; likewise should stay bounded 
both as as C — > —00 and as ( — > 0. 

We are now in a position to analyze the behavior of the 
adjoint eigenmodes; we will report the analysis in some 
detail, as the various elements form important ingredi- 
ents of the derivation of a moving boundary approxima- 
tion in the next section. C* simplifies again considerably 



in the positive £ region due to the fact that bo vanishes 
identically there, so it is again of advantage to split the 
region of integration into two, £ < with C* given by 
(60), and C > for which 



C* = 



-Vo-^+no 



a 2 



-n Q 

- v 



For ( > $2 has to solve the homogeneous ODE: 



<9C 2 



• v - 



0. 



(64) 



(65) 



This equation has two independent solutions, a constant 
and an exponential that diverges for increasing £. Hence 
the boundary condition that "J 2 remains bounded imme- 
diate gives the solution 



*2 = ipo = constant, £ > 0. 



(66) 



Moreover, the differential equations for \& 2 implied by the 
zero-eigenvalue equation 



C* 



*2 



= 



(67) 



shows that ^2 has to be continuous and have a contin- 
uous derivative at £ — 0. Hence, when we construct the 
eigenmodes on the left half space £ < 0, the \& 2 compo- 
nent has to obey the boundary condition 



*2(C = o) = Va), 9* 2 /aci c =o = 0. 



(68) 



Since b vanishes identically for £ > 0, we need to 
know ^1 only in the region £ < 0. As we stated above, 
because the functions b\ that we consider all vanish as 
£ j 0, the definition of the adjoint operator does not 
imply a boundary condition on ^\{C = 0) as long as it 
does not diverge. A straightforward analytical investiga- 
tion of the equation near £ — shows that in general 
will, with a finite slope, approach a finite value as £ | 0, 
and that in general it has a higher order singular term 

_ (_£)(l + D/fc)_ 

We now turn to the behavior as £ — > —00. In this limit, 
n — > and b — > 1, so C* reduced from (60) to 




a? 







(69) 



It is easy to verify that as £ — ► —00, there are three 
possible types of non-divergent zero mode solutions, 



a voC/D 



$0) 



~ e 



(70) 



Here A± = (t>o±y / t> M^t)/2, so that the mode 
deed converges towards the left; the other mode allowed 
by the linear equations, e x -'>, on the other hand, diverges 
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towards the left, and hence is forbidden by the boundary 
conditions. 

The mode "F^ is very special — it is immediately ver- 
ified from (60) that 



ipo 
V'o 



= for all ( (tpo const.), 



(71) 



not just for £ — > — oo. In other words, the constant mode 
"F^ 1 ) is an exact adjoint zero mode for all ( < 0. 

If we integrate or forward towards increas- 
ing (, each trajectory in the phase space of the ODE 
is uniquely determined (apart from an overal amplitude, 
as the equations are linear). Hence, if we follow either 
\|/(2) or ,jf(3) towards £ = 0, the derivatives |f =0 
and c^'F^ 3 -' |^ = o wm i n general be nonzero. As we have 
however seen above, in order that the full eigenmodes on 
the whole real axis remain bounded also for ( — > oo, the 
^-component needs to have a zero derivative at £ = 
— see Eq. (68). Each separate eigenmode does not obey 
this requirement, but by the linearity of the equation it 
will always be possible to construct one unique linear 
combination <F( 2 ) of "F^ and which does have zero 
derivative. This solution constitutes the second adjoint 
zero eigenmode of the problem. Like the trivial eigen- 
mode "F^, it can be extended smoothly to the required 



(2) 



const, behavior for £ > 0. 
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FIG. 12. Zero right eigenvector ^' 2 ' of the adjoint oper- 
ator C* for k = 2 and D — 0.3. The 6-like component 
approaches a finite value at £ = with a finite slope; it has 
generally a higher order singular term ~ (_^)( 1 + D / fe ). 

Fig. 12 shows the adjoint zero eigenmode 'F^ 2 ' of £* 
for k = 2 and D = 0.3, obtained numerically from the 
ODE's with a shooting method. The qualitative behavior 
of the components is in agreement with the above discus- 
sion, and is independent of the values of the parameters. 

- (2) 

Note that the 6-like component ' approaches a finite 
value at C — with a finite slope; this behavior is also 
found for arbitrary parameters, while as is easily verified 



there generally is a parameter-dependent subdominant 
singular term proportional to \(\( 1+D / k ) . 

To obtain the onset of lateral instability the adjoint 
zero eigenmode has to be convoluted with the transla- 
tional mode (d^bo,d^no) according to (62). Which of 
the two zero modes should we use? The trivial ad- 
joint zero mode "F^ expresses change of velocity under 
reparametrization, and is related to conservation law in 
our system [see the discussion after Eq. (26)]; this will 
become more clear in the next section. It does not play a 
role for the onset of instability; wc will therefore ignore it 
here [33] and use "F^ to evaluate (62). A change of sign 
of the numerator, which marks the onset of instability, 
is obtained for k = 0.5,1,2 and 3 and shown in Fig. 3 
as diamonds. In the figure, these values are also com- 
pared to the D c determined by the numerical dispersion 
relation shown in Fig. 3 as filled dots. The agreement is 
very good, as it should; we have also checked that away 
from this line, a fit of the small-g behavior of the dis- 
persion relation leads to values of the second derivative 
of lo at q = which are consistent with the solvability 
formula. These results thus confirm the consistency of 
our full stability calculation and the solvability expres- 
sion for the critical line in the phase diagram and the 
small q behavior of the growth rate uu. 



IV. SHARP INTERFACE FORMULATION 

A moving boundary approximation or sharp interface 
formulation is appropriate when the width of the front 
or interface is much smaller than the typical length scale 
of the pattern and when the dynamics of the pattern oc- 
curs through the motion of these interfaces. The moving- 
boundary approximation amounts then to treating these 
fronts as mathematically sharp interfaces or boundaries 
by taking their width to zero and integrating out their 
internal degrees of freedom. There are three important 
assumptions underlying such an approximation, namely 
(a) that there is a separation of length scales, (b) that 
there is a separation of time scales between the motion 
of the front as a whole and its internal dynamics, and (c) 
that the internal dynamics of the front is determined by 
the nonlinear front region itself, so that the solvability- 
type integrals are dominated by the contributions from 
the finite region, and hence do not diverge. The latter 
condition is violated in practice only for special types 
of fronts propagating into a linearly unstable state, so- 
called "pulled" fronts [34,35]; our fronts are not of this 
type (they are of the "pushed" type, in this terminol- 
ogy), so we focus our analysis on the length and time 
scale requirements (a) and (b). 

As we saw in section III, the planar front width is fi- 
nite, even for D — > at fixed k. Moreover, even though 
for k — > 00 the 6-field rises over an exponentially small 
distance behind the singularity line, both the b and n field 
even then only approach their asymptotic values over a 
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distance of order unity. In this sense, even in this limit 
the front width W remains finite. Of course, we can al- 
ways choose to investigate fronts whose curvature k is 
small in the sense that kW <C 1. For these, a mov- 
ing boundary approximation should be accurate; we do 
find indeed below that the sharp interface formulation 
we derive for the present problem is consistent with the 
result of the dispersion relation of section III for q small 
enough that qW <C 1. However, whether such a moving 
boundary approximation applies at all dynamically rele- 
vant length scales, is another matter! We already know 
from the analysis of the dispersion relation in section III 
that in the left part of the D-fc-phase diagram modes 
up to q = q c are unstable, and that q c is generally fi- 
nite, except close to the critical line D = D c (k). Hence, 
q c W remains finite as well, and so there is no obvious 
limit where a moving boundary approximation becomes 
asymptotically correct on all dynamically relevant length 
scales. Nevertheless, we find that in practice the sharp 
interface formulation that we develop is rather accurate 
in a significant portion of the phase diagram. Since the 
present problem has some unusual and new aspects, we 
focus again on the essential structure and intuitive argu- 
ments, rather than mathematical rigor. 



A. Sharp Interface Formulation of the Problem 

The simplest case to consider to guide our intuition 
is the limit D <C 1. As we discussed in connection 
with Fig. 5, in this regime the bacterial density field 
approaches its asymptotic value on the finite scale W, 
while the n-diffusion field in front of the bacterial front 
decays on a length scales l n = 1/vq, which diverges as 
D — > since vq ~ D. A sharp interface formulation is 
then based on the idea that we view the bacterial front on 
the "outer" scale £ out on which the patterns and diffusion 
fields vary in the presence of the moving boudary, which 
is treated as a sharp line of zero thickness. The dynamics 
of the fields on the "inner" scale [36,37] of the front (W), 
and the way in which this dynamics responds to changes 
in the fields on both sides of this inner zone, is trans- 
lated into boundary conditions for the outer fields in the 
interfacial formulation. Formally, the moving boundary 
consists of taking the limit 5 = W/£ ou t — ► 0. 

Normally, a sharp interface formulation or moving 
boundary approximation is based on applying the the- 
ory of matched asymptotic expansions [28,37]. Here, the 
situation is somewhat unusual: on the right side of the 
inner region (the interfacial transition zone where essen- 
tially the nutrient is consumed by the bacteria) we al- 
ready have a sharply defined boundary, the singular line 
where b vanishes. At this side, we therefore do not have 
a matching problem, instead we already have two bound- 
ary conditions for the n field, namely the requirements 
that the value of n and the gradient of n are continous 
at the singular line where b vanishes — this follows di- 



rectly from the dynamical equation (7), since the "reac- 
tion term" — nb is continuous. On the left side of the 
interfacial zone, on the other hand, the b field varies con- 
tinuously and we do have a true matching problem. 

In our present case, the "outer field" on the front side 
of the interface is simply the n field in the whole region 
to the right of the singular line where b vanishes; hence 
there the nutrient field n obeys a simple linear diffusion 
equation: 

r 6=o, 

front side "outer" eqs.: < (72) 
l = V 2 n. 

It is useful to introduce a suitable curvilinear coordinate 
system in which £ = coincides with the singular line 
where b vanishes. In the sharp interface limit, the line 
£ = then also coincides with the position of the mov- 
ing interface. We furthermore identify the region ahead 
of the front as the + side of the interface where £ > 0, 
and use a superscript + to indicate values of the outer 
field extrapolated from the outer region towards the line 
£ = 0: n + = lini£ m n(£), Vn+ = lim^o Vn(£). As we 
already mentioned, n and its gradient should be continu- 
ous at this line, and we can therefore write the boundary 
conditions as 

limra(£) = n+, limVn(f) = Vn+. (73) 

We now turn to the behavior on the back side of the 
front, the — side. We have seen that in the bacterial 
front, the nutrition field decays exponentially fast to zero 
towards the left on a length scale of order unity; hence, in 
the — outer region behind this interfacial zone, we have 
n w for the nutrition field. The bacterial field is close 
to 1 there. Therefore, we take the dynamics of the 6-field 
into account there, by writing b = 1 + Ab and linearizing 
in the outer field Ab, 

( §*> = DV 2 Ab, 
back side "outer" eqs.: < (74) 
[ n = 0. 

Note that the outer equations (72) and (74) have been 
written in the laboratory frame, not in a co-moving 
frame, since in the case of nontrivial patterns, there is 
no single relevant co-moving frame. 

What are the boundary conditions on the — side of 
the front? According to the matching prescription, the 
inner field expanded in the outer variables on the back 
side should equal the outer field expanded in the inner 
variables [37]. Extrapolating the inner field in the outer 
variables towards the — side means that we investigate 
the 6-profile to the left towards — oo. As Figs. 5 and 6 
illustrate, on the inner scale W the b field rapidly ap- 
proaches a constant value. Although these figures are 
made for planar fronts, the analysis below shows that 
this continues to hold for weakly curved fronts and that 
the appropriate matching conditions are 
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lim^-oo 6(0 = 1 + A6-, 
matching condition: ^ (75) 
limf—oo V6(£) = 0. 

Here Ab~ is the value of the outer field Ab extrapolated 
from the outer — region towards the interface. Note that 
the second condition that the gradient vanishes, is also 
consistent with the matching prescription: if we assume 
that the outer field Ab varies on the outer scale X = 5r 
with 5 = W I l ou t then the outer gradient of Ab rewritten 
in terms of the inner variable vanishes in the limit 5^0. 

Now that we know how to connect the inner fields to 
the outer ones — on the left side of the inner region 
through the matching conditions (75), on the right side 
through boundary conditions (73) — we are ready to de- 
rive the effective boundary conditions in a sharp-interface 
formulation. One easily convinces oneself that in order 
to get a well-posed moving boundary problem with the 
above outer dynamical equations and matching condi- 
tions, one needs effectively three boundary conditions re- 



lating A& 



and Vn + and the interface velocity and 



curvature. To derive them, we imagine that the front is 
weakly curved with curvature k such that kW oc 6 -C 1. 
Since we identified the line £ = with the line where b 
vanishes, £ is a local co-moving coordinate with speed v 
in the direction perpendicular to the front. In this weakly 
curved system, we then have on the inner scale 



db _ D 

dt ~ ' ~ 

dn 
~dt 



Q2 h k+1 

k+i &e + ' u 

dn 



-Dnb k 



db 



- nb, 



d 2 n , , dn 



(76) 
(77) 



Following standard practice, we now ignore the time 
derivatives on the left (taken in the co- moving frame). 
This amounts to an adiabaticity assumption, the assump- 
tion that the change in the pattern and hence in the front 
speed and profile, take place on time scales much longer 
than the relaxation time of the front (assumption (b) dis- 
cussed in section IV. A above) . Technically, it means that 
the solution stays always close to a uniformly translating 
solution in the curved coordinate system, and our goal 
now is to calculate the changes in the velocity perturba- 
tively. Indeed, let us write v = Vo + v\, where v a is the 
velocity of the planer solution and v\ the change in veloc- 
ity due to the curvature and the fact that the outer field 
n is changed slightly from the planar solution; similarly 
we write b = b a + b[ and n = n + n[ so that b[ and n[ 
are the deviations in the fields from the planer solutions 
(we use the prime to emphasize the difference from the 
perturbations used in the linear stability analysis). Upon 
linearization about the planar front solution, we then get 
the equation 



C 



- Vl -Dnb^ 

—V\ — K 




(78) 



where C is the same operator introduced before in (38), 
written now in terms of the variable £. This equation 



again calls for applying the solvability condition. We 
have already seen in section III.C that the operator C 
has a number of adjoint zero cigenmodes. There is a 
subtle difference between the present analysis and the 
one of section III.C, however, which is crucial for ob- 
taining the proper boundary conditions. In the stability 
analysis of section III.C we worked with functions de- 
fined on the whole interval (—00,00); this implied that 
the mode constructed for negative £ needed to have zero 
derivative of the ^2 component at £ = 0, and this re- 
duced the number of proper adjoint zero eigenmodes to 
two. Here, however, we are doing perturbation analy- 
sis only in the inner region, which in the inner variable 
is the semi-infinite interval (— 00, 0]. We therefore do 
not need to require now that the adjoint zero mode has 
d^2/d^ = o — 0, and hence we now have three rather 
than two admissable adjoint zero modes! These lead to 
the three necessary equations that become the boundary 
conditions sought for in the sharp interface formulation. 
Moreover, because we now work on a semi-infinite inter- 
val, we get boundary terms at £ = from the partial 
integrations necessary to have the operator work on the 
adjoint functions [39]: 



J — ( 



d£ (*i,* 2 ) -C 



-f 

J — ( 



d£ 



*2 



-^i(-oo)v Ab - 



8^2 



9£ 



(n+-n+) 



+ * 2 (0) [V 5 n+ - V 5 n+ + v (n+ - n+)] . (79) 

Here we have used the boundary and matching condi- 
tions (75) and (73) for the deviations b[ and n[ from the 
planar values bo and no- Note that there are no bound- 
ary terms in the field b[ at £ = 0, since these are all 
proportional to 6§, and 6§(£ — ► 0) vanishes according to 
(23). 

The three boundary conditions now straightforwardly 
follow by taking the left inner product of the three zero 
modes with equation (78) together with (79). The be- 
havior of the three adjoint zero modes of C* on the left 
for £ < has been discussed already in section III.C. The 



first one is simply = ty^ 1 ' = constant. In this case 
we immediately obtain 



(i) 



-v Ab + {V 5 n + - V e n+) + v a (r 



-LA 



rf£( fa+DKb^-^ + fa+K]- 



dnQ 

3£ ' L ^ ' " J &f 
= «i[l-nJ] + /s[D/(fc +!)-<]. (80) 



We note that this equation is essentially a type of con- 
servation equation in a weakly curved frame — indeed, it 
can also be obtained by an analysis similar to the deriva- 
tion of the conservation equation (26) by adding the two 
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equations (76), and (77) and integrating, ignoring the 
temporal derivates for a quasi-stationary front solution 
in the co-moving frame. This is the reason for our earlier 
remark in section III.C that the constant left zero mode 
^f^ 1 ) is related to conservation. 

The second and third boundary conditions are ob- 
tained from the two other adjoint zero modes , j( 2 ) and 
, j( 3 ) of C* , discussed in section III.C; the general form is 
similar to the one shown for a ^ 2 - ) for a particular choice 
of parameters in Fig. 12, except that the derivative of the 
n-like component does not vanish at £ = 0. These zero 
eigenmodcs of C* can only be evaluated numerically, but 
the form of the boundary condition is simply the same 
for both: with i = 2, 3 we get 

*«(0)((Vin+-V ? n+) 




In the sharp interface interpretation, equations (80) and 
(81) are interpreted as the boundary conditions that re- 
late the change in the local field n[ and its gradient at 
the interface (relative to those for a planar moving front) 
to the local change in velocity of the interface and the 
local curvature. For a general pattern, the derivative of 
■n' x with respect to £ on the left of these equations has 
to be interpreted as derivative in the direction normal to 
the interface [40] . As we discussed above, these equations 
are precisely the three boundary conditions necessary to 
get, together with the outer equation (72), a well-posed 
moving boundary problem. 



that the gradient of the bacterial field b does not ap- 
pear, makes these bacterial growth fronts most like the 
so-called "one-sided crystal growth models" , describing 
situations where the diffusion on the back side is absent 
(e.g., in "directional solidification", the diffusion of im- 
purities in the solid on the back side of the interface is 
usually negligible in comparison with the diffusion in the 
liquid on the front side [7]). 

We can make this observation more precise as follows. 
Note that the boundary condition (80), the one that ex- 
presses conservation, is the only one which involves Ab~ . 
Hence we can solve the full dynamical problem by work- 
ing only with the outer n-field and the two boundary con- 
ditions (81) — in other words, the outer n-field together 
with these boundary conditions constitute a closed prob- 
lem that is sufficient to describe the dynamics of the mov- 
ing interface. Once this is determined, one can use the 
conservation condition (80) to determine Ab~ and from 
there analyze the behavior of the 6-field on the back with 
the outer equation (74): in the sharp interface limit the 
b-field becomes completely slaved to the interface motion! 

We already anticipated in section III.B that the cur- 
vature correction term (the effective surface-tension-like 
term) should be of order D for small D. This is fully 
confirmed by our analysis: all the curvature terms in the 
boundary conditions (80) and (81) cither explicitly in- 
volve a term D, or a term proportional to no which, as 
we saw in section II, is proportional to D too for small 
D. 



C. Applicability to the various regimes 



B. Interpretation of the sharp interface formulation 

It is useful to pause for a moment to reflect on the 
structure of the boundary conditions. First of all, note 
that they all involve terms proportional to the gradient 
of the nutrient diffusion field, and not to the bacterial 
density field. The presence of these gradients of the nu- 
trient field on the front side of the interface could have 
been expected from the numerical observations that the 
bacterial growth fronts arc unstable for small enough D. 
It is well know [7] that such interfacial instability arise 
for diffusion limited growth problems where the growth 
velocity of the interface is proportional to the gradient 
of the driving field that "feeds" the interface. The fact 



The simplest way to test the accuracy of a moving 
boundary formulation is by comparing the dispersion re- 
lation from the moving boundary problem with the dis- 
persion relation obtained from the full model as discussed 
in section III. A. The two outer equations are linear dif- 
fusion equations of the standard form, while the bound- 
ary conditions are (by construction) also linear. Conse- 
quenty, the stability of the planar solution of the sharp 
interface problem follows the standard stability problem 
as discussed in [7] in which the growth or decay of small 
single-mode perturbations about the planar interface so- 
lution is determined. We will therefore not report it ex- 
plicitly here. 
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FIG. 13. Comparison between dispersion relation obtained 
through a sharp interface approach (full line) and through 
direct numerical linear stability analysis (dots), for the case 
k = 2, D = 0.001. Note that for q < 1 the results from the 
moving boundary approximation agree very well with those 
of the full linear stability calculation, as it should. The reason 
for the difference between the two curves for q values of order 
unity is discussed in the text. 

In Fig. 13 we illustrate such a comparison for a typical 
case in the small D regime, where the moving boundary 
approximation is expected to work best since the diffu- 
sion length in the nutrient + region. The figure, which 
is for the case k = 2 and D = 0.001 confirms that for 
small q the dispersion relation of the moving boundary 
approximation (full curve) essentially lies on top of the 
one derived from the full problem (symbols). This is as 
it should, since for small q clearly qW <C 1, so that the 
condition for the moving boundary approximation to be 
accurate is fullfillcd. The overal shape of the dispersion 
relation of the moving boundary approximation is actu- 
ally quite close to the exact one, but for larger q there 
clearly are some quantitative differences, even for this 
small value of D. This discrepancy is in our view due 
to what we discussed before, the fact that the range of 
unstable wavenumbers for this case is finite (q c ~ 0.8), 
while the interface width W is finite too, so that q c W 
does not approach zero as D — > 0. 

Even though the moving boundary approximation 
therefore does not become formally correct in this limit 
(in the sense that the correction terms can not be made 
arbitrarily small by taking D sufficiently small) , it clearly 
does quite well in practice for these parameter values. 
Probably this is due to the fact that W is still relatively 
small compared to the wavelength A c = 2n/q c corre- 
sponding with the marginal wavenumber q c : If we take 
W ~ 2 we get W/X c « 1/4, so even though we can not 
make this ratio arbitrarily small by sending D — > 0, it ap- 
pears to be small enough in practice to make the sharp 
interface formulation work well. What may also play a 
role is that for problems with nonlinear diffusion like this 
one, the response of the interfacial zone to perturbations 
is mostly determined by the rather thin zone where b is 
small; we have not attempted to substantiate this intu- 
itive idea, however. 



A similar observation holds for the timescalcs. As 
Fig. 9 illustrates, the maximum growth rate Lu m of the 
most unstable mode is proportional D and hence vq for 
small D; the proportionality uj m ~ v q m is also con- 
sistent with the Mullins-Sekerka dispersion relation (57) 
discussed in section III. A. The internal relaxation time 
t front of the front itself is expected to be of order W/v, 
hence Lo m T ~ q m W remains finite in the limit D — > 0: 
there is no full separation of timescales either. 

Also for k S> 1 and D of order unity, the present ap- 
proximation works generally very well, since on the one 
hand the diffusion length in the nutrient zone ahead of 
the front is large (as vo ~ 1/fc for k large), while on 
the other hand the interfacial zone tends to becomes rel- 
atively small, even though it does not appear to go to 
zero — see Fig. 6 and the discussion at the end of section 
II. We have also investigated the possibility whether in 
the limit k 3> 1 another approximation might be possi- 
ble, one in which there are three zones, an outer region 
in front of the interface where b = as we had above, 
a very thin zone (of exponentially small width, see sec- 
tion II) where b quickly rises to a value close to 1 while 
n hardly changes, and a region behind this zone where 
Ab = b — 1 is already small and where n decays to zero. 
We have not been able to make this approximation work 
to our satisfaction, basically because we have not been 
able to match the thin zone properly to the region behind 
it. 

For values D and k of order unity, but not too close to 
the stability boundary D c (k), the discrepancy between 
dispersion relation of the moving boundary approxima- 
tion that we have derived above and the exact dispersion 
relation is bigger than in Fig. 13 for small D. This is to be 
expected, since in this limit diffusion length in the outer 
n-region is of the same order as the interface width, and, 
moreover, the diffusion in the bacterial region is more 
important. Nevertheless, the order of magnitude of the 
growth rate and the range of unstable wavenumbers is 
right. 

Finally, we note that since a long-wavelength instabil- 
ity occurs upon decreasing D below D c {k), we expect 
that just to the left of this line, the dynamics can be de- 
scribed by the so-called Kuramoto-Sivashinsky equation 
[41,42]. In fact, since the line D c (k) is very straight for 
k larger than 1, it may be well be the problem also sim- 
plifies in the limit D^oo,k^oo,D/k fixed. We have 
not attempted to study this limit or to give an explicit 
derivation of the Kuramoto-Sivashinksky equation near 
the boundary. 



V. SUMMARY AND OUTLOOK 

We have shown that a nonlinear diffusion coefficient 
and a simple bilinear autocatalysis is sufficient to gen- 
erate a long wave length instability as long as the dif- 
fusion constant obeys D < D c , where D c depends on 
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the nonlinearity. We hope that these results will be of 
help in sorting out to what extent the present class of 
models does describe the real bacterial growth problems. 
To do so, one would of course have to be able to map 
the bacterial growth properties onto the effective diffu- 
sion coefficient in this model. If this can be done, the 
most clear test with the aid of the present results would 
be to see whether the interfacial instability becomes sup- 
pressed once the effective bacterial density becomes too 
large. 

It would also be of interest to extend numerical sim- 
ulations like those of Golding et al. [10] and those of 
Kitsunczaki [22], whose parameter values are indicated 
by crosses in Fig. 3: As we discussed in the introduction, 
these numerical results appear to contradict our analyti- 
cal results, but this may be due to finite size effects and 
the limited time of the simulation. 

In addition to providing a starting point for further 
studies of these models for bacterial growth, we have been 
able to develop a new type of linear stability calculation 
which applies to the general class of fronts with singular 
behavior of the fields as a result of nonlinear diffusion. 
Our methods therefore opens up the possibility to study 
other systems as well, like the vortex fronts [19,20] we 
mentioned in the introduction. In addition we were able 
to reformulate the reaction diffusion problem with non- 
linear diffusion away from the stability boundary D c (k) 
in the form of a free boundary type problem; our analysis 
shows that in the present model bacterial growth fronts 
are closest to those described by the so-called one-sided 
model in the theory of crystal growth [7]. 

Of course, our results are far from the final answer 
on these type bacterial growth models: the knowledge 
that the planar front is unstable is only the first (though 
crucial) step towards understanding the actual evolving 
patterns, which are determined by nonlinear effects. In 
addition, within the context of understanding the bacte- 
rial growth problem, the question remains to what extent 
models with nonlinear diffusion suffice to capture the im- 
portant growth dynamics. 

Clearly, we have studied only the simplest variant of 
such models, by leaving out the death term, which ap- 
pears to be important for the morphology [22,2], as well 
as lots of effects which are important for a more realistic 
model for bacterial colony growth, like the sporulation 
of bacteria already mentioned in the introduction: Com- 
puter simulations have shown that in order that branches 
can form, the sporulation term — \xb has to be present in 
(6). 

We have not tried to study how the critical line D c (k) 
approaches the origin as k I 0; this is an interesting tech- 
nical question, but one which probably is of limited rel- 
evance for understanding the bacterial growth problem. 
In fact, the behavior near the origin in the D-fc-phase di- 
agram is very singular and hence sensitive to changes in 
the model: the nonlinear diffusion as well as finite cutoff 
effects as well as changeing the bilinear reaction term bn 
to Ifn change the mathematical behavior dramatically. 



Finally, we want to draw attention to an open mathe- 
matical question — at least for us. In a solvability type 
analysis, the boundary conditions one normally imposes 
on the adjoint fields follow from the requirement that 
(*(£$)) = ((£**)$) for all dynamically relevant func- 
tions <3? [see Eq. (63)]. However, in the derivation of our 
moving boundary approximation, we have operated dif- 
ferently: instead of imposing boundary conditions on the 
adjoint functions, we have written out the terms from the 
partial differential equations explicitly, and used the left 
zero modes on the half space (-co, 0] that we already 
knew to obtain the sought for boundary conditions for 
the physical fields ni! Clearly, the equations obtained 
this way follow necessarily from the original differential 
equations in the weakly curved frame, but this line of 
reasoning is mathematically different in spirit from the 
usual Fredholm alternative (solvability theory). We do 
not know — nor could we find — the mathematical the- 
ory behind this approach which appears to be new and 
very powerful for problems with a singular line. 
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